Equilibrium Bundle Size of Rodlike Polyelectrolytes with Counterion-Induced 

Attractive Interactions 
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Multivalent counterions can induce an effective attraction between like-charged rodlike polyelec- 
trolytes, leading to the formation of polelectrolyte bundles. In this paper, we calculate the equilib- 
rium bundle size using a simple model in which the attraction between polyelectrolytes (assumed to 
be pairwise additive) is treated phenomenologically. If the counterions are point-like, they almost 
completely neutralize the charge of the bundle, and the equilibrium bundle size diverges. When 
the counterions are large, however, steric and short-range electrostatic interactions prevent charge 
neutralization of the bundle, thus forcing the equilibrium bundle size to be finite. We also consider 
the possibility that increasing the number of nearest neighbors for each rod in the bundle frustrates 
the attractive interaction between the rods. Such a frustration leads to the formation of finite size 
bundles as well, even when the counterions are small. 

PACS numbers: 61.20.Qg, 61.25.Hq, 82.35.Rs, 82.70-y 
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The mean-field Poisson-Boltzmann (PB) theory of 
electrostatic interactions predicts that two identical 
macromolecules in any salt solution will repel each other 
1.]. However, the presence of multivalent counterions 
can actually induce an attraction between like-charged 
polyelectrolytes (PEs). This has been experimentally 
observed for several different PEs, including double- 
stranded DNA HQ, F-actin 00, microtubules 0,11], 
and the fd, M13, and tobacco mosaic viruses 0,0- Com- 
puter simulations of both homogeneously c harg ed rods 
E0E3EI1 and realistic DNA molecules [H III un- 
ambiguously show that attractive interactions can arise 
solely from counterion correlations not included in PB 
theory. Several theories that take these correlations into 
account - including perturbative expansions of PB the- 
ory 0,^3; structural-correlation theory [IHIl7lll^ |. and 
strong-coupling theory [13 _ obtain an attractive inter- 
action between two rods, ft is still a matter of discussion, 
however, as to which of these theories is the most appro- 
priate description of the correlation-induced attraction 
seen in experiments and simulations. Furthermore, it is 
unknown whether the interactions between multiple rods 
is pairwise additive or not [17], l20i |2l| . 

Under experimental conditions in which the interaction 
between PEs is attractive, the PEs typically form dense, 
ordered bundles of a well-defined size, rather than pre- 
cipitating into a PE-rich phase 0,0, H, H • In this 
paper, we theoretically investigate the thermodynamic 
stability of these bundles (if bundle growth is not limited 
thcrmod yna mically, then it must be limited by kinetic 
barriers |2l], |22j, |23j ) . We assume that the attractive in- 
teractions are pairwise additive, but do not specify the 
precise nature of the counterion correlations. Rather, 



we simply introduce a phenomenological parameter 7 to 
characterize the attractive energy between two PEs in a 
bundle. 

Consider, then, an aqueous solution of volume V with 
N identical rodlike PEs of length L, radius ao, and a 
uniform linear charge density — eAo (the aggregation of 
flexible PEs has been considered in 22] ). We treat the 
aqueous solution as a uniform dielectric with dielectric 
constant e, and ignore all image-charge effects (i.e. we 
assume that the dielectric constant of the rods is f» e). 
Positive monovalent and g-valent counterions, as well as 
negative monovalent co-ions, are present; the entire sys- 
tem is charge neutral and in chemical equilibrium with 
a salt bath. It is well known that a highly charged PE 
will have a larg e number of counterions "condensed" near 
its surface . A condensed multivalent ion can become 
strongly correlated both with ions condensed on the same 
rod and with ions condensed on nearby rods. The latter 
correlations give rise to the effective attraction between 
rods, which in turn leads to the formation of PE bundles 
of some size at equilibrium. For simplicity, we assume 
that only multivalent ions can enter inside the bundle 
(the effects of competitive binding with monovalent ions 
will be discussed in a future paper |23|V Also, we assume 
that the solution of rods is dilute, i.e. the volume frac- 
tion cf> = Nna^L/V -C 1. As a result, we can employ the 
cell model, where each bundle and its surrounding ions 
are enclosed in a Wigner-Seitz (WS) cell, and interac- 
tions between cells are ignored. We work in the long-rod 
limit (L — > 00), so that the translational entropy of the 
bundles is negligible. Finally, we assume that the equi- 
librium distribution of bundle sizes is sharply peaked, so 
that all bundles in the system have approximately the 
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same size. Given these assumptions, the free energy can 
be written as a sum of four terms: 



N 

PF=— [pF ent 



• PF ES + f3F attr ] + Nf3F c , 



(1) 



where M is the number of rods in each bundle and 
f3 = l/k B T, k B being Boltzmann's constant and T the 
temperature. (3F ent includes the entropy and the chemi- 
cal potential terms for all of the ions in a single WS cell; 
(3Fes is the total mean-field electrostatic energy for a WS 
cell; [3F at tr is the total correlation-induced attractive en- 
ergy for a single bundle; and (3F corr is the correlation 
energy for ions condensed on a single rod. 

To calculate (3F ent , we assume that the ions inside the 
bundle are uniformly distributed throughout the volume 
available inside the bundle. The ions outside the bundle, 
on the other hand, are treated with a modified Debye- 
Huckel (DH) theory similar to Manning's counterion con- 
densation theory |2J]. It is well known that DH theory 
breaks down near highly charged surfaces |25|: in partic- 
ular, it does not properly account for counterion conden- 
sation. To correct for this, we allow ions to condense into 
a Stern layer surrounding the bundle. We assume that 
the ions inside the layer are uniformly distributed. The 
ions outside the layer are described by the distribution 
functions n s (x), where s = ±1, q is the ionic species. For- 
mally, Debye-Huckel theory can be derived by expand- 
ing the free energy to quadratic order in the difference 
n s (x) — n s 25], where n s is the bulk concentration of ion 
species s. Assuming the ions are point-like and discard- 
ing terms that contribute constants to the free energy, 
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where X q L is the number of multivalent ions condensed 
on each rod in the bundle, a b ML = ttL (a 2 — a 2 ,) is the 
volume available to the ions inside the bundle, XfL is the 
number of s-valent ions in the Stern layer, and a s tL w 
2ttRwL is the volume of the layer (w <C R is the width 
of the layer). 

In order to calculate the mean-field electrostatic energy 
of the system, we model each bundle as a homogeneously 
charged cylinder of radius R « a\M (so that its vol- 
ume is approximately equal to the bundle volume), and 
the surrounding Stern layer as a uniform surface charge 
distribution (i.e. we set w = 0). That is, — en b (x) = 
—9(R — r)eX/Tra 2 and en s t(x) = 9{R — r)eX s t/2irR are 
the charge distributions of a bundle and its Stern layer, 
respectively, where — eX = —e (Ao — qX q ) is the renormal- 
ized linear charge density of one rod in the bundle and 
eX st = e (Af* + qX s q l — Al') is the linear charge density of 
the Stern layer. The electrostatic free energy is given by 
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where l B = e 2 /ek B T is the Bjerrum length (l B ~ 7.lA 
in water) and entot(x) = e J2 S sn s {x)6(r — R) + en st (x) — 
en b (x) is the total charge distribution for a WS cell. 

The correlation-induced attractive interactions lead to 
the formation of a "bond" of energy E bon d — -~lk B TL 
between each pair of neighboring rods in the bundle. If 
we assume that the rods in the bundle are packed in a 
hexagonal array, as has been experimentally observed 0, 
0- S- 01 - then, to a very good approximation for all M > 2, 

PF attr = /3E bond B = - 7 L (3M - 3.6VM) , where B is 
the number of bonds in a bundle. The first term is the 
bulk attractive energy; the second term is an effective 
surface tension that arises from the fact that the rods 
on the bundle surface have less neighbors than the bulk 
rods. Thus, we can see that each bundle in our model is 
equivalent to a homogeneously charged cylinder with an 
effective surface tension in the presence of counterions. 
This is very similar to the Rayleigh instability j2(| of a 
charged water droplet in the presence of counterions (2lJ • 
The number of ions inside each bundle and its sur- 
rounding Stern layer, as well as the distribution of ions 
outside of these regions, minimize the total free en- 
ergy (3F (subject to the constraint of overall charge 
neutrality). For the ion distributions n s (x), this mini- 
mization yields the expected DH distributions, n s (x) = 
n s [1 — sijj(x)], where the dimensionless electrostatic po- 
tential ip(x) is given by 
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Here, —eXtot = —eXM + eX s t is the total linear 
charge density of the bundle and its Stern layer, k 2 — 
47tZb Yls s2n si an d K u (x) is the modified Bessel func- 
tion of the second kind of order v. Clearly, both 7 and 
0F corr should depend on X q . However, since the pre- 
cise nature of the counterion correlations is not speci- 
fied in our model, we ignore this dependence; that is, we 
calculate X q at the mean-field level. Discarding all con- 
stant terms, the free energy eq. Q can be written as 
j3F = N [a 2 Ti/R 2 + T 2 L] where T Y is given by the final 
three terms of eq. © and 

3.67a l B X 2 R 2 l B a 2 X 2 ot K ( K R) 
•Fi = — ^ 1 — 1 tv»^ ; — ( 5 ) 
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The minimization conditions for X q and A** are, respec- 
tively: 
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FIG. 1: Free energy 0F(R) for a = Iran, a = 
lAnm, w — Iran, A^ 1 = 0.17 nm, n\ = 10 mM, q = 3, 7 = 
1.4nm _1 , and n 5 = 0.1 /jM (solid line), 10 )j,M (dashed line), 
and ImM (dotted line). Inset: Xtot(R) for the above param- 
eter values with n q = 10 fj.M. The dashed line indicates the 
asymptotic value of Xtot (see text). 



FIG. 2: Free energy f3F(R) [2jj for finite size ions inside the 
bundle, with 7 = 1.7nm _1 and A* = O.OlAo (the remai ning 
parameter values are the same as in Fig. |Q. Inset: f3F(R) |2fl| 
for frustrated attractive interactions, with £ = 0.1, <f)max = 2, 
and 7 = 1.7 nm -1 (the remaining parameter values are the 
same as in Fig. 
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If we solve eqs. © and (0), we can see that A tot ap- 
proaches a constant at large i? (see inset of Fig. [[). 
Indeed, in the limit il — * 00, the second term in eq. 
© - which is due to the electrostatic self-energy of the 

A ° ^ and 



A 



bundle - dominates, causing A — > , 2a p2 In 

tot ~^ 1+kw ■ Thus, as the bundle grows, additional 
counterions condense inside the bundle, so that the total 
charge of the bundle remains constant. 

The equilibrium bundle radius R eq is given by the value 
of R that minimizes the free energy. Fig. ^ shows the 
free energy as a function of R for various values of n q . 
We can see that as soon as the attractive energy is strong 
enough to induce bundle formation, R, 



eq 



00. This is 

a result of the large number of counterions that enter 
each bundle, which causes the entropic and electrostatic 
resistance to bundle growth to be weak: in the limit R — > 

00, l/M(/3F E s + 0F ent ) 1/R 2 . This resistance is 

overwhelmed by the attractive energy (3F attr /M ~ 1/R, 
thus causing the equilibrium bundle size to diverge. 

As stated above, these results do not take into account 
the dependence of 7 and l3F corr on X q . It is important 
to note, however, that the asymptotic results |A| ~ 1/R 2 
and (3F ~ 1/R — 1/R 2 hold for any such dependence 
(even those which cause overcharging): when the bundle 
is large, the long-range mean-field self-energy of the bun- 
dle dominates over the short-range correlation energy, 
causing the renormalized charge density of the bundle to 
be small and the resistance to bundle growth to be weak. 

In the model discussed above, the density of ions inside 
the bundle can in principle be arbitrarily high. In reality, 
however, both steric and short-range electrostatic inter- 
actions limit the ion density inside the bundle. Steric 



interactions prevent the ion density inside the bundle 
from exceeding the close packing density. In addition, 
our model underestimates the strong repulsion between 
the ions inside the bundle when their density is high, 
because the charge of these ions is smeared out over the 
entire bundle in calculating (3Fes ■ This additional short- 
range electrostatic repulsion effectively increases the size 
of the ions in the bundle. If we treat the ions inside the 
bundle as finite size particles with an effective volume 
v eff , then (3F ent -> 0F ent +ML(A-X q ) ln(l-A 9 /A)+A g , 
where A » cvfc /v e f f is the maximum number of ions per 
unit length that can condense on a single rod. This adds 
a term - ln(l - X q /A) to the RHS of eq. © which di- 
verges as X q — > A, thus forcing X q < A for any bundle 
size R. If A* = Ao — qA > 0, then the asymptotic result 
A ~ 1/R 2 no longer holds; rather, A — > A* for large R. 
In order for A* > 0, the counterions must exceed a min- 
imum effective size set by the spacing between the rods 
in the bundle (we assume that this spacing is set by the 
minimum of the short-range correlation-induced attrac- 
tion energy, which is independent of the bundle size). For 
DNA, where 1/A = ljA , a = 10A, and a ps 14A in 
the presence of trivalent cobalt hexammine 131, this 
implies the effective radius of the ions 6 e ff > 7A. When 
A* > 0, the self-energy of the bundle - in particular, the 
second term in eq. (jSJ - diverges at large R, as shown in 
Fig. |3 A local minimum in (3F(R) is obtained only if the 
surface tension dominates the free energy for some values 
of R. As shown in Fig. |2 this is not the case at small 
multivalent salt concentrations, as the total attraction 
is not large enough to overcome the electrostatic repul- 
sion between rods. As n q increases, however, more ions 
enter inside the bundle and reduce the electrostatic re- 
pulsion for small bundles, thus creating a local minimum 
in (3F{R) that is primarily determined by the balance of 

the first two terms in eq. R cq 
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tice that the bundling transitions is discontinuous, and 
that the bundle size should be independent of further 
increases in n q , as has been recently observed for micro- 
tubule bundles |2s|. 

Up to this point, we have assumed that the attraction 
between neighboring rods is independent of the bundle 
size. There are several mechanisms, however, that can 
frustrate the bonds between rods as the bundle grows. 
For short-ranged, pairwise additive interactions, only a 
fraction of ions in a narrow "contact stripe" between 
the two rods become correlated with one another [l7| . 
When rod-rod dimers form, the size of the contact stripe 
is maximized. If the size of this stripe is large enough, 
then the bond energy will decrease as the bundle grows 
and the rods have to divide the ions in their condensed 
layers equally among all of their bonds. Alternatively, 
a non-uniform PE charge distribution results in a rela- 
tive orientation that minimizes the electrostatic repul- 
sion between two neighboring rods. Achieving the op- 
timum orientation between a rod and all of its near- 
est neighbors may cost energy or be physically impossi- 
ble; in either case, the bond energy effectively decreases 
as the bundle grows. Indeed, it has been experimen- 
tally observed that F-actin filaments undergo twist dis- 
tortions when forming bundles to reduce the electro- 
static repulsion between neighboring rods e 5] . If we write 
PF a ttr = -7-LM^(^), where B/M = b is the aver- 
age number of bonds per rod (b < 3 for hexagonally 
ordered bundles), then cj>(b) — b for unfrustrated inter- 
actions. To encapsulate the effects of bond frustration, 
we chose <f)(b) to be a hyperbola that approaches the 
asymptote <p = b for small b and <f> = 4> m ax for large b, 

<f>(b) = i (6+ <f>rnax) - § \<f>max ~ &| yj 1 + JpZ^^f ■ In 

other words, the total attractive energy gained for each 
rod in the bundle saturates as b increases; by adjusting 
4>max and £, we can control the saturating value and rate 
of saturation, respectively. As shown in the inset of Fig. 
121 a local minimum in (3F(R) can be obtained by choos- 
ing certain values of (frmax and £. Unlike the minimum 
obtained with finite-size ions, this effect cannot lead to 
arbitrarily large bundles; rather, the onset of frustration 
must occur at a sufficiently small bundle size, or the en- 
tropic and electrostatic resistance to bundle growth will 
not be large enough to prevent infinite bundles. 

In summary, we have introduced a mean-field model to 
calculate the equilibrium bundle size of highly charged, 
rodlike polyelectrolytes in the presence of multivalent 
counterions. For point-like counterions and a short-range 
attraction that is independent of the bundle size, many 
counterions enter inside the bundle, causing the short- 
range attraction to overwhelm the resistance to bundle 
growth and leading to infinite bundles at equilibrium. In 
order to obtain a finite equilibrium bundle size, then, ei- 
ther the electrostatic self-energy must be enhanced or the 
attractive energy suppressed for large bundles. The for- 



mer can be accomplished if the short-range interactions 
between the ions inside the bundle are large enough to 
prevent the neutralization of large bundles by counterion 
condensation; the latter, if the interactions between rods 
in the bundle become frustrated as the bundle grows. 
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